# -*- coding: utf-8 -*-

"""Created on Sat Apr 11 17:39:16 2020 @author: L JL
   zq 优化于 2021.10.31

"""
import numpy as np
import pandas as pd
from common.config import torqueConfig as Config



'''
# 输入：(x,y)：某节点坐标，(X,Y)：N个数据点坐标，N：数据点数量, d:截断半径
# 输出：从小到大排列的距离，shape:[N, 2],其中[:,0]:从小到大的距离，[:,1]相应点index
'''
def D_radius(x, y, X, Y, N):
    ind_mat = np.zeros((N, 2))  # ind_mat[权重,数据点号]
    for i in range(N): #400个数据点用时 0.009s
        s = ((x - X[i]) ** 2 + (y - Y[i]) ** 2) ** 0.5
        ind_mat[i][0] = s
        ind_mat[i][1] = i
    # 按距中心点距离s排序
    ind_mat = ind_mat[np.lexsort(ind_mat[:, ::-1].T)]  #::-1表示翻转该维度的数据，相当于左右列对掉; np.lexsort(次排序数组,主排序数组)得到排序索引，400个数据点用时0.001s
    return ind_mat


'''
#计算权重
#输入：（截断半径，节点坐标x0,y0，数据坐标x[N],y[N]）
#输出：权重，标量
'''
def W_mat(d, x0, y0, x, y):
    s = (((x - x0) ** 2 + (y - y0) ** 2) ** 0.5) / d
    W = ((2/3)-4*s**2+4*s**3) * (s<=0.5) + ((4/3)-4*s+4*s**2-(4/3)*s**3) * (0.5<s) * (s<=1)
    return W


'''
#
#输入：W:np.array[N],N个数据点对该node节点的权重，P:np.matrix[N,3]:3个基在N个数据点出的值，M:基数量，N：数据点数量
#输出：主方程式中的A矩阵
'''
def A_mat(W, P, M):  # 系数矩阵
    W = np.array([W for _ in range(M)])
    A = np.multiply(P.T, W)@P
    return A

'''
#
#输入：u:np.array[N]:N个数据点的Z坐标，W:np.array[N]:N个数据点对该node节点的权重，P:np.matrix[N,3]:3个基在N个数据点出的值，M:基数量，N：数据点数量
#输出：np.array[M]主方程中额B矩阵
'''
def B_mat(u, W, P):
    B = np.array(np.multiply(W, u)@P).reshape(-1)
    return B



class SurfaceDecorator:

    def fit(fun):
        print("surfaceDecorator.fit 进行中 ..........")
        def wrapper(points):
            x, y, z = points[Config.WindSpeed].values, points[Config.RotateSpeed].values, points[
                Config.ActPower].values
            #print("x1: ", x)
            #print("y1: ", y)
            #print("z1: ", z)

            # np.array形式的测量点坐标，行向量
            X = np.array(x)
            Y = np.array(y)
            Z = np.array(z)
            Xmax = (np.max(X))
            Xmin = (np.min(X))
            Ymax = (np.max(Y))
            Ymin = (np.min(Y))
            N = len(X)  # 总点数
            M = 3  # 基底矢量数

            '''
            主方程式：A*a = B
            [(p1,p1)...(p1,pj)...(p1,pm)] [a1_node]   [(u_p,p1)]
            [......                     ] [.......]   [.......]
            [(pi,p1)...(pi,pj)...(pi,pm)]*[aj_node] = [(u_p,pi)]
            [......                     ] [.......]   [.......]
            [(pm,p1)...(pm,pj)...(pm,pm)] [am_node]   [(u_p,pm)]
            其中，pi表示第i个基；ai_node表示pi系数；u_p表示点p测量值；变量p表示基，下标p表示测量点。
                  (pi,pj)=sum_所有N个测量点( w(x_node-x_p)*pi(x_p)*pj(x_p) )
                  (u_p,pi)=sum_所有N个测量点( w(x_node-x_p)*u_p*pi(x_p) )
            '''
            P = np.mat(np.zeros((N, M)))  # 该node节点处的M个基在N个测量点处的值
            u = np.array(Z)  # 测量点的Z坐标，行向量
            # W，A，B，a，d_mat 变量在后面都会作为左值出现，因此写在这里是多余的，但是有助于理解其数据结构，因此保留下来。
            W = np.mat(np.zeros((N, 1)))  # N个测量点在该node节点处的权重w？？
            A = np.mat(np.zeros((M, M)))  # 上面的P矩阵
            B = np.mat(np.zeros((M, N)))  # ？？？
            a = np.mat(np.zeros((M, 1)))  # 上面的系数矩阵，aj_node列矩阵，好像没有用上！！！！！！！！！
            d_mat = np.mat(np.zeros((N, 1)))  # dataZ = []？？，好像没用上！！！！！！！！
            dataX = np.arange(Xmin, Xmax, Config.deltaX)  # X方向格点坐标(风速）
            dataY = np.arange(Ymin, Ymax, Config.deltaY)  # Y方向格点坐标（风向）
            #print("Xmin, Xmax: ", Xmin, Xmax)
            #print("Ymin, Ymax:", Ymin, Ymax)

            #优化2：由于P[N,M]的值是跟所选的node节点位置无关的，因此，将其赋值从下面循环中提取出来
            for n in range(0, N):
                P[n, 0] = 1
                P[n, 1] = X[n]
                P[n, 2] = Y[n]  # 选用的是一次曲面拟合

            points_fited = pd.DataFrame(columns=[Config.WindSpeed, Config.RotateSpeed, Config.ActPower])
            for i in dataX:
                for j in dataY:
                    d = Config.neighborRadiu  # 权重w的截断半径
                    ind_mat = D_radius(i, j, X, Y, N)
                    # print(ind_mat)
                    if ind_mat[3, 0] <= d:  # 如果以该node节点为中心，半径d内部至少有3个点，进行曲面拟合，计算节点Z坐标
                        try:
                            W = W_mat(d, i, j, X, Y)
                            A = A_mat(W, P, M)
                            B = B_mat(u, W, P)

                            c = np.linalg.solve(A, B)
                            # dataZ.append(c[0]+c[1]*i+c[2]*j)
                            dataZ = c[0] + c[1] * i + c[2] * j  # 该node节点处的拟合值
                            #print("A", end=' ')

                        except:  # 如果拟合不成功，取最近的4个数据点Z坐标平均值
                            ind_Zsum = 0
                            for ind in range(0, 4):
                                ind_Zsum += Z[int(ind_mat[ind, 1])]  # dataZ.append(ind_Zsum/4)
                            dataZ = ind_Zsum / 4
                            print('B', end=' ')
                                # print(dataZ)

                    else:  # 如果以该node节点为中心，半径d内部 不足 3个点，取最近的4个数据点Z坐标平均值
                        ind_Zsum = 0
                        for ind in range(0, 4):
                            ind_Zsum += Z[int(ind_mat[ind, 1])]
                            # dataZ.append(ind_Zsum/4)
                        dataZ = ind_Zsum / 4
                        print('C', end=' ')
                            # print(dataZ)

                    if dataZ < 0:  dataZ = 0  # 如果格点附近没有测试点，可能出现该点拟合值比较奇异，这里限制一下
                    points_fited = points_fited.append(pd.DataFrame(
                        {Config.WindSpeed: [i], Config.RotateSpeed: [j], Config.ActPower: [dataZ]}),
                                                       ignore_index=True)
                    #print(i, j, dataZ)


            #points_fited.to_csv(Config.MLSSurfaceFitProcess, index=False)

            #print("points_fited: ", points_fited)

            return fun(points_fited)
        return wrapper

